library(ggplot2)
library(patchwork)
'%!in%' <- function(x,y)!('%in%'(x,y))

# set to desired working directory
# setwd("")

# Start log file for script -----------------------------------------------
TeachingDemos::txtStart("./log files/figureA23.txt")

# Load Data ---------------------------------------------------------
# individual data
load("./data/impdat20.Race.WeekWeight.rdata")
load("./data/impdat20.Age.WeekWeight.rdata")
load("./data/impdat20.Arab.WeekWeight.rdata")
load("./data/impdat20.Asian.WeekWeight.rdata")
load("./data/impdat20.Disab.WeekWeight.rdata")
load("./data/impdat20.GendCar.WeekWeight.rdata")
load("./data/impdat20.GendSci.WeekWeight.rdata")
load("./data/impdat20.NatAm.WeekWeight.rdata")
load("./data/impdat20.Pres.WeekWeight.rdata")
load("./data/impdat20.Religion.WeekWeight.rdata")
load("./data/impdat20.Sexuality.WeekWeight.rdata")
load("./data/impdat20.Weapons.WeekWeight.rdata")
load("./data/impdat20.Weight.WeekWeight.rdata")

# * Age ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Young_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Age.WeekWeight, 
                    subset = race5 == "White" & date_mdy != "2020-05-25",
                    weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Age"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,54,107)] <- "2020-01-01"

# Plotting
age_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Arab ------------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.OtherPeople_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Arab.WeekWeight, 
                    subset = race5 == "White" & date_mdy != "2020-05-25", 
                    weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Arab"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
arab_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Asian -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Asian.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Asian"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
asian_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Disability --------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Abled_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Disab.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Disability"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
disability_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# * Gender-Career ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Male_Career_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendCar.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Career"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
gendcar_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Gender-Science ----------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Male_Science_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.GendSci.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Gender-Science"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
gendsci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Native American ---------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.White_American_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.NatAm.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Native American"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
natam_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# President ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Trump_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Pres.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "President"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
pres_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Religion ----------------------------------------------------------------

# ^ Christianity-Islam ----------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Christianity_Ci_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Islam"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
relig_ci_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# ^ Christianity-Judaism --------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Christianity_CJ_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Christianity-Judaism"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
relig_cj_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())



# ^ Judaism-Islam ---------------------------------------------------------

# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Judaism_Ji_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Religion.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Judaism-Islam"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
relig_ji_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Sexuality ---------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Straight_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Sexuality.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Sexuality"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
sexuality_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Weapons -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Black_Weapons_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weapons.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weapons"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
weapons_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())




# Weight -----------------------------------------------------------------
# ** Model ----------------------------------------------------------------
# GF is week 21 (a Monday)
m1_ideo_dummy <- lm(D_biep.Thin_Good_all ~ 
                      as.factor(week)*ideo_mod + 
                      as.factor(week)*ideo_con + 
                      age_cat4 + edu_cat4 + sex + 
                      census_region + broughtwebsite2,
                    data = impdat20.Weight.WeekWeight, subset = race5 == "White" & date_mdy != "2020-05-25", weights = WeekWeights)
summary(m1_ideo_dummy)

# ** Figure ----------------------------------------------------------
title <- "Weight"
weeks <- names(table(m1_ideo_dummy$model$`as.factor(week)`))

pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")
pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

preds_lib <- predict(m1_ideo_dummy, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo_dummy, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo_dummy, pdat_con, se.fit = T)

# combine into DF to plot
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01"

# Plotting
weight_plot <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2) +
  geom_linerange(aes(ymin = pred - 1.96*se, 
                     ymax = pred + 1.96*se),
                 lwd = .75) +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              # method = lm,
              aes(x = week, y = pred)) +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "", title = title) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 12),
        legend.text = element_text(size = 12),
        strip.text = element_text(size = 14),
        strip.background = element_blank())


# Race IAT for Placebo Comparison -----------------------------------------
m1_ideo <- lm(D_biep.White_Good_all ~ 
                as.factor(week)*ideo_mod + 
                as.factor(week)*ideo_con + 
                age_cat4 + edu_cat4  + sex + 
                census_region + broughtwebsite2,
              data = impdat20.WeekWeight, weights = WeekWeights, 
              subset = race5 == "White" & date_mdy != "2020-05-25")

weeks <- names(table(m1_ideo$model$`as.factor(week)`))

# Data Frames for Predictions
pdat_lib <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

pdat_mod <- data.frame(week = weeks,
                       ideo_mod = 1,
                       ideo_con = 0,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

pdat_con <- data.frame(week = weeks,
                       ideo_mod = 0,
                       ideo_con = 1,
                       race5 = "White",
                       age_cat4 = "18-29",
                       edu_cat4 = "Post-Grad",
                       sex = "f",
                       census_region = "South",
                       broughtwebsite2 = "Assignment school/work")

weeks <- names(table(m1_ideo$model$`as.factor(week)`))
preds_lib <- predict(m1_ideo, pdat_lib, se.fit = T)
preds_mod <- predict(m1_ideo, pdat_mod, se.fit = T)
preds_con <- predict(m1_ideo, pdat_con, se.fit = T)

# combine into DF
p_dat <- data.frame(week = weeks,
                    ideo = c(rep("Liberal", length(weeks)),
                             rep("Neutral", length(weeks)),
                             rep("Conservative", length(weeks))),
                    pred = c(preds_lib$fit, preds_mod$fit, preds_con$fit),
                    se = c(preds_lib$se.fit, preds_mod$se.fit, preds_con$se.fit))
p_dat_fill <- data.frame(week = unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)],
                         ideo = c(rep("Liberal", length(unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)])),
                                  rep("Neutral", length(unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)])),
                                  rep("Conservative", length(unique(impdat20.WeekWeight$week)[which(unique(impdat20.WeekWeight$week) %!in% weeks)]))),
                         pred = NA,
                         se = NA)
p_dat <- rbind(p_dat, p_dat_fill)
p_dat$ideo <- factor(p_dat$ideo,
                     levels = c("Liberal",
                                "Neutral",
                                "Conservative"))
p_dat$week <- as.Date(paste(2020, p_dat$week, 1, sep="-"), "%Y-%U-%u")
p_dat$week[c(1,length(weeks)+1, length(weeks)*2+1)] <- "2020-01-01" # from 00 weeks


# Plot
race_fig <- ggplot(p_dat, aes(x = week, y = pred)) +
  geom_point(size = 2, color = "grey") +
  geom_linerange(aes(ymin = pred - 1.65*se, 
                     ymax = pred + 1.65*se),
                 lwd = .75, color = "grey") +
  geom_smooth(data = subset(p_dat, week < "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  geom_smooth(data = subset(p_dat, week >= "2020-5-25"), 
              span = .8,
              aes(x = week, y = pred), se = FALSE, color = "black") +
  facet_grid(ideo~., scales = "free_y") +
  theme_bw() +
  labs(y = "Predicted D-Score",
       x = "") +
  scale_y_continuous(labels = numform::ff_num(zero = 0, digits = 2)) +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 20),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 20),
        strip.background = element_blank())

lhs_fit <- ggplot_build(race_fig)$data[[3]]
lhs_fit$y[which(lhs_fit$PANEL == 1)]
lhs_fit$y[which(lhs_fit$PANEL == 2)]
lhs_fit$y[which(lhs_fit$PANEL == 3)]

rhs_fit <- ggplot_build(race_fig)$data[[4]]
rhs_fit$y[which(rhs_fit$PANEL == 1)]
rhs_fit$y[which(rhs_fit$PANEL == 2)]
rhs_fit$y[which(rhs_fit$PANEL == 3)]


race_preds1 <- data.frame(iat = "race",
                          out = "loess",
                          Liberal = rhs_fit$y[which(rhs_fit$PANEL == 1)][1] - lhs_fit$y[which(lhs_fit$PANEL == 1)][length(lhs_fit$y[which(lhs_fit$PANEL == 1)])],
                          Neutral = rhs_fit$y[which(rhs_fit$PANEL == 2)][1] - lhs_fit$y[which(lhs_fit$PANEL == 2)][length(lhs_fit$y[which(lhs_fit$PANEL == 2)])],
                          Conservative = rhs_fit$y[which(rhs_fit$PANEL == 3)][1] - lhs_fit$y[which(lhs_fit$PANEL == 3)][length(lhs_fit$y[which(lhs_fit$PANEL == 3)])])
race_preds2 <- data.frame(iat = "race",
                          out = "predicted",
                          Liberal = p_dat[which(p_dat$week == "2020-05-25" & p_dat$ideo == "Liberal"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" & p_dat$ideo == "Liberal"), "pred"],
                          Neutral = p_dat[which(p_dat$week == "2020-05-25" & p_dat$ideo == "Neutral"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" & p_dat$ideo == "Neutral"), "pred"],
                          Conservative = p_dat[which(p_dat$week == "2020-05-25" & p_dat$ideo == "Conservative"), "pred"] - p_dat[which(p_dat$week == "2020-05-18" & p_dat$ideo == "Conservative"), "pred"])
race_preds <- rbind(race_preds1, race_preds2)
race_preds <- reshape::melt(race_preds, id = c("iat", "out"), variable_name = "ideo")


# Plot -----------------------------------------------------
weeks <- unique(p_dat$week)

obs_diffs <- as.data.frame(array(NA, c(14, 4)))
names(obs_diffs) <- c("iat", "Liberal", "Neutral", "Conservative")

placeb_diffs <- as.data.frame(array(NA, c(1, 5)))
names(placeb_diffs) <- c("iat", "week", "Liberal", "Neutral", "Conservative")

x_diffs <- as.data.frame(array(NA, c(length(weeks), 5)))
names(x_diffs) <- c("iat", "week", "Liberal", "Neutral", "Conservative")

plot_names <- grep("_plot", ls(), value = T)
plot_names <- plot_names[which(plot_names != "facet_plot")]

# loop across plot objects to pull predicted D-scores for observed differences and pseudo differences
for(i in 1:length(plot_names)){
  # metadata storing
  x_diffs$iat[1] <- gsub("_plot", "", plot_names[i])
  x_diffs$week[1] <- weeks[1]
  
  # pull observed George Floyd Difference
  obs_diffs$iat[i] <- gsub("_plot", "", plot_names[i])
  
  lhs_fit <- ggplot_build(get(plot_names[i]))$data[[1]]
  rhs_fit <- ggplot_build(get(plot_names[i]))$data[[2]]
  
  obs_diffs$Liberal[i] <- rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == "2020-05-18")[1]]
  obs_diffs$Neutral[i] <- rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == "2020-05-18")[1]]
  obs_diffs$Conservative[i] <- rhs_fit$y[which(rhs_fit$PANEL == 3)][which(p_dat$week == "2020-05-25")[1]] - lhs_fit$y[which(lhs_fit$PANEL == 3)][which(p_dat$week == "2020-05-18")[1]]
  
  # other differences but for the treatment period
  for(j in 2:length(weeks)){
    x_diffs$iat[j] <- gsub("_plot", "", plot_names[i])
    x_diffs$week[j] <- weeks[j]
    
    if(weeks[j] <= "2020-05-18"){
      x_diffs$Liberal[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$Neutral[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$Conservative[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 3)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 3)][which(p_dat$week == weeks[j-1])[1]]
    }
    if(weeks[j] == "2020-05-25" | weeks[j] == "2020-06-01"){
      # empty placeholder
    }
    if(weeks[j] > "2020-06-01"){
      x_diffs$Liberal[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 1)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 1)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$Neutral[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 2)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 2)][which(p_dat$week == weeks[j-1])[1]]
      x_diffs$Conservative[j] <-  rhs_fit$y[which(rhs_fit$PANEL == 3)][which(p_dat$week == weeks[j])[1]] - lhs_fit$y[which(lhs_fit$PANEL == 3)][which(p_dat$week == weeks[j-1])[1]]
    }
  }
  
  # store results
  placeb_diffs <- rbind(placeb_diffs, x_diffs)
}

# Reformatting and cleaning up
placeb_diffs <- placeb_diffs[-1,]
placeb_diffs <- reshape::melt(placeb_diffs, id = c("iat", "week"), variable_name = "ideo")
obs_diffs <- reshape::melt(obs_diffs, id = c("iat"), variable_name = "ideo")

# Getting p-values 

# Data frame to store p-values
p_values <- obs_diffs

for(i in 1:nrow(p_values)){
  # store test statistic
  ob_val <- obs_diffs$value[i]
  # extract appropriate simulations for reference distribution
  ref_dist <- subset(placeb_diffs, iat == p_values$iat[i] & ideo == p_values$ideo[i],
                     select = value)
  ref_dist <- ref_dist$value[which(!is.na(ref_dist$value))]
  # store number of simulations
  n_sim <- length(ref_dist)
  # two-tailed p-value
  p_values$value[i] <- sum(abs(ref_dist) > abs(ob_val))/n_sim
}
p_values$x <- -.15
p_values$y <- 15
p_values$p <- paste0("p=", round(p_values$value, 2))

unique(placeb_diffs$iat)

labs <- c("Age", "Arab", "Asian", "Disability",
          "Gender-\nCareer", "Gender-\nScience",
          "Native\nAmerican", "President", 
          "Religion\nChristian-\nIslam", "Religion\nChristian-\nJudaism",
          "Religion\nJudaism-\nIslam", "Sexuality", "Weapons", "Weight")

placeb_diffs$iat <- factor(placeb_diffs$iat, labels = labs)
p_values$iat <- factor(p_values$iat, labels = labs)
obs_diffs$iat <- factor(obs_diffs$iat, labels = labs)


# Plotting
ggplot(placeb_diffs, aes(x = value)) +
  geom_histogram(fill = "grey") +
  geom_vline(data = obs_diffs, aes(xintercept = value)) +
  geom_text(data = p_values, aes(label = p, x = x, y = y)) +
  facet_grid(iat~ideo, scales = "free_x") +
  scale_x_continuous(labels = numform::ff_num(zero = 0, digits = 2), limits = c(-.35, .35)) +
  theme_bw() +
  labs(x = "Predicted IAT Score Difference", y = "Count") +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 16),
        axis.title = element_text(size = 14),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        strip.text = element_text(size = 16),
        strip.background = element_blank())
ggsave("./figures/FigureA23.pdf", height = 16, width = 8)

# End log file for script -------------------------------------------------
TeachingDemos::txtStop()